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Abstract 

We propose a practical method to solve the random-phase approximation (RPA) in the self- 
consistent Hartree-Fock (HF) and density-functional theory. The method is based on numerical 
evaluation of the residual interactions utilizing finite amplitude of single-particle wave functions. 
The method only requires calculations of the single-particle Hamiltonian constructed with inde- 
pendent bra and ket states. Using the present method, the RPA calculation becomes possible with 
a little extension of a numerical code of the static HF calculation. We demonstrate usefulness and 
accuracy of the present method performing test calculations for isoscalar responses in deformed 
20Ne. 
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I. INTRODUCTION 



The mean-field theory with a density-dependent effective interaction has been an essential 
tool to understand nuclei. Thanks to the high performance computing, it is now becoming 
the most promising tool for quantitative description of nuclear structure in medium-to-heavy 



nuclei tU, y] • The nuclear self-consistent mean-field theories are analogous to to the density- 
functional theory in condensed matter. A current major goal is constructing a universal 
energy-density functional, which is able to describe ground and excited states in nuclei and 
nuclear matter. This is also urgently needed for predicting and interpreting new data from 
the next generation of radioactive beam facilities. 

In order to describe dynamical properties in nuclear response to external fields, the 
random-phase approximation (RPA) is a leading theory applicable to both low-lying states 
and giant resonances Sj. The RPA is a microscopic theory which can be obtained by 
linearizing the time-dependent Hartree-Fock (TDHF) equation, or equivalently, the time- 
dependent Kohn-Sham equation in the density-functional theory. The linearization produces 
a self-consistent residual interaction, v = 6'^E[p\/6p'^, where E and p are the energy-density 
functional and the one-body density, respectively (Sec. [TTl)- The standard solution of the 
RPA is based on the matrix formulation of the RPA equation, which involves a large num- 
ber of particle- hole matrix elements of the residual interaction, Vph',hp' and Vpp'^hh'- Since the 
realistic nuclear energy functional is rather complicated, it is very tedious and difficult to 
calculate all the necessary matrix elements. It is, therefore, the purpose of the present paper 
to present an alternative method of solving the RPA equations, in which we deal with only 
the single-particle Hamiltonian, h[p]. 

Although there are numerous works on the HF-plus-RPA calculations, because of the 
complexity of the residual interactions, it has been common in practice to neglect some 
parts of the residual interactions. The RPA calculations with full self-consistency are be- 
coming a current trend in nuclear structure studies, however, they are essentially only for 



spherical nuclei at present 



The applications to deformed nuclei are very few, but 



have been done___for the _Skyrme energy functional using the three-dimensional mesh-space 



representation 



111 ]. See Sec. I in Ref. j4| for a current status of these studies. 



The basic idea of the present method is ana 



a time-dependent manner (real-time method) lOl . 



ogous to linear-response calculations in 
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12 1 . In the real-time method, the 
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time evolution of a TDHF state involves only the action of the HF Hamiltonian, h[p{t)], 
onto single-particle orbitals, \ipi{t)) {i = 1, - ■ ■ ,A). Although the real-time method is very 
efficient for obtaining nuclear response in a wide energy range, its numerical instability 



caused by zero modes was a problem for the linear-response calculations Zero-energy 
modes related to symmetry breaking in the HF state are easily excited, which often prevents 
the calculation of the time evolution for a long period. Therefore, it is desirable to develop 
a corresponding method in the frequency (energy) representation. 

This paper is organized as follows. A new approach to the solution of the linear response 
equation, "finite amplitude method", is presented in Sec. [Tll In Sec. IIIH using the Bonche- 
Koonin-Negele (BKN) interaction [13], we check the accuracy of solutions obtained with 
the present method. We also investigate the zero-energy components in calculated strength 
functions. Then, the conclusion is summarized in Sec. IIVI 



II. LINEAR RESPONSE THEORY 



The RPA equation is known to be equivalent to the TDHF equation in the small- 
amplitude limit. We recapitulate how the RPA equation is derived from the small-amplitude 
TDHF equation, that will help explanation of our basic idea. 



A. TDHF and linear response equation 

The HF Hamiltonian, h[p] = SE[p\/Sp, is a functional of one-body density matrix, p, 
which satisfies the condition, = p, that the state is expressed by a single Slater determi- 
nant. The stationary condition is 

[h[p],p]=0, (1) 

which defines the HF ground state density p = Pq. Hereafter, the static HF Hamiltonian 
is simply denoted as Hq = /i[po]) ^^nd h = 1 is used. When a time-dependent external 
perturbation is present, the time evolution of the density, p{t), follows the TDHF equation, 

if^p{t)=[h[p{t)]+VUt), Pit)]. (2) 

Using this p{t), the expectation value of a one-body operator F is obtained as (F) = 
tr{Fp(t)}. Provided that the perturbation is weak, we may linearize Eq. ((21) with respect 
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to Vcxt{t) and Sp(t) defined by 

pit)=po + 5p{t). (3) 
This leads to a time-dependent linear-response equation with an external field, 

'jt^P(t) = N' ^P(^)] + [^ext(t) + Sh{t),po], (4) 

where 6h{t) is a residual field induced by density fiuctuations, 

(5) 



Shit) ^ — ■ 5p{t) = J2 



P=PO 

It should be noted that 5h{t) has a linear dependence on Sp(t). As we will see in Eq. ( fl^ . if 
we adopt the natural basis diagonalizing ho, the summation can be restricted to the particle- 
hole {p > A, u < A) and hole-particle {p > A, u < A) components. Now, we decompose the 
time-dependent Sp{t) into those with fixed frequencies: 

Spit) = E iv^Pi^)^"''' + V*Sp\iv)e"-'}. (6) 

UJ 

The external and induced fields are also expressed in the same way. 

Sh{t) = J2{vSh{uj)e-'^' + ri*5h\uj)e'^'} , (7) 

UJ 

KxtW = Yl {vVU^)e-^^' + v*vU^)e'^'} ■ (8) 

UJ 

Here, we have introduced a small dimensionless parameter rj. 6h{uj) may be written as 
6h{uj) = 6h/6p ■ 6p{uj). Note that the transition density, the external field, and the in- 
duced field in the ^-representation, Sp^u), Vcxti^^), and 6h{Lj), are not necessarily hermitian. 
Substituting these into the linearized TDHF equation (jlj), we obtain the linear- response 
equation in the frequency representation, 

u 5p{uj) = [ho, 5p{uj)] + [Vext(t^) + 5h{uj),po\. (9) 

This is the equation we want to solve in this paper. 

When the frequency u is equal to an RPA eigenfrequency Un, there is a non-zero solution, 
6pn, of Eq. Qj with V^xt = 0. These are called normal modes and are orthogonal to each 
other. The normalization is given by 

H[^pLSPn']Po} = {<^o\[Spi,Spn']\'^o) = Snn', (lO) 



where |$o) indicates the HF ground state. From Eq. (fTOj) . it is obvious that, in order to 
normahze the transition density 5p„, it must be non-hermitian. When uj = Un, the nucleus 
is truly excited by Vext(i^), and we cannot determine the magnitude of Sp{u!n) because Sp(t) 
increases in time. If 5p{uJn) is a solution of Eq. ([H]), 5p{uJn) +c6pn with an arbitrary constant 
c is a solution too. 

So far, the linear-response equation has been expressed in terms of the one-body density 
operators. The density-matrix formulation is simple and easy to manipulate, however, in 
practical calculations, it is convenient to introduce single-particle (Kohn-Sham) orbitals. For 
systems with A particles, the TDHF describes the one-body density using A single-particle 
orbitals, \i/ji(t)), 

A A 

Pit) = po = El'^^)^'^^!- (11) 

1=1 1=1 

It is an advantage of the TDHF, that the time evolution is described by occupied orbitals 
only, {l^/'i)} with i = 1, ■ ■ ■ , A. The static orbitals are normally chosen as eigenstates of the 
HF Hamiltonian, 

ho\^,) = e,\^^), (12) 

which can be divided into two categories; occupied (hole) orbitals, {(pi} (i = 1, ■ ■ ■ ,A), for 
which we use indexes and unoccupied (particle) orbitals, {0m} (m = yl + 1, ■ ■ ■ ), for 

which we use indexes m,n, - ■ ■ . In the linear approximation, we have 

^Pit) = J2{\^^){mt)\ + \5Ut)){^^\}, (13) 
i 

where \ipi(t)) = (|(/)j) + |(5?/'j(t)))e~"'* and it is linearized with respect to \SilJi(t)). The 
condition, p(t)^ = pit), leads to 

Spij = 6pmn = 0, i,j<A, m,n>A, (14) 
+ (5^,10,) = 0. (15) 

The second equation is nothing but the orthonormalization condition for single-particle 
orbitals, {\ipi{t))} (i = !,■■■ ,A). 

Transforming 6p{t) into Sp{uj) in Eq. ([6]), we must make ket and bra states independent, 
because Sp^u) is not hermitian. This is related to the fact that the RPA equation is described 
by forward and backward amplitudes, X{uj) and Y{uj). 

5p{cu) = J2{\M^)){<P^\ + \<P^){Y^i^)\}■ W 
i 
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This is equivalent to the Fourier decomposition of the time-dependent single-particle orbitals, 

= E {v\X^{u))e~'^' + r^*\Y,{u))e^-'}. (17) 

Since only the particle-hole matrix elements of Sp{uj) are non-zero, seen in Eq. fll4l) . we can 
assume that the amplitudes, \Xi{u!)) and \Yi{u!)), can be expanded in the particle orbitals 
only; 

\X,iuj)) = E |0^^)X™(^), \Y,iuj)) = \<t>n.)Y:,,iuj). (18) 

m>A m>A 

If we take particle-hole and hole-particle matrix elements of Eq. (Q, we can derive the 
well-known RPA equation in the matrix form; 





Here, the matrices, A and B, and the vectors, / and g, are defined by 

dh 



(19) 



I 



dp, 



nj 



^.), (20) 



Bmi,nj — (0m I 



P=PO 

dh 



■jn 



(21) 



P=PO 



fmi{uj) = (0m|V;xt(t^)|0i), 9mii^^) = (0i | Kxt (t^) 1 0m) • (22) 

For the normal modes that are homogeneous solutions of Eq. f|T9l) . the orthogonality and 
normalization are expressed as 

Ef x-(n') _ ^^{n)*^^{n') \ _ r fr)o\ 

y^mi ^mi ^ mi ^ mi j ~ ^nn' ^ \^'^) 



m,i 



which is equivalent to Eq. flTU]) . This is a standard matrix formulation of the RPA equation. 



In practical applications, the most tedious part is ca^ 



culation of matrix elements of the 



residual interactions in Ami,nj and Bmi,nj- In Ref. jl4|, a numerical method to solve the 
RPA equation in the coordinate space is proposed, and the similar approaches are used in 
realistic applications using the Skyrme interaction . In those works, one does not need 
to calculate the particle orbitals, however, the residual interaction must be evaluated in 
the coordinate-space representation. In Sec. Ill B\ we propose an alternative, even simpler 
approach to a solution of the linear-response equation The method does not require 
explicit evaluation of the residual interaction, 6h/6p. 



B. Finite amplitude method 

Multiplying both sides of Eq. ([9]) with a ket of hole states we have 

uj\X,{uo)) = {ho - e,)\X,{uj)) + Q {KxtM + Sh{uj)} (24) 

where Q is a projection operator onto the particle space, Q = 1 — Another 
equation can be derived by multiplying a bra state with Eq. (j9]). 

uj{Y,{u)\ = -{Y,{u)\{ho - e.) - {<Pi\ {Kxtl^} + Sh{uj)} Q. (25) 

These are formally equivalent to the RPA equation in the matrix form of Eq. ([T^ . 

The essential idea of our new numerical approach is as follow: Equations (1241) and (125]) 
require operations of the HF Hamiltonian in the ground state, ho, and the induced fields, 
5h{uj) and 6h'^{u!). Since ho is obtained by the static HF calculation, a new ingredient for 
the RPA calculation is the latter two. The conventional approach is to expand 6h{uj) in the 
linear order as Eq. ([5]), then to solve the RPA equation in a matrix form. In this paper, 
instead of performing the explicit expansion, we resort to the numerical linearization. Now, 
let us explain how to achieve it. 

The time-dependent self-consistent Hamiltonian, h{t), is a functional of one-body density 
that is represented by occupied A single-particle states; h[p{t)] = hlipit)]. In the linear 
approximation, 

h[po + 6p{t)] = h[(f) + 6ij{t)] = ho + 6h{t), (26) 

the induced field can be written as 6h(t) = h[(f) + Sip(t)] — ho- In the frequency representa- 
tion, the story becomes slightly more complicated, because 5h{uo) and 5h'^{(jj), are no longer 
hermitian. In this case, we should regard h[p] as a functional of 2 A single-particle states 
(independent bra and ket), and I'j/'j), i = 1, ■ ■ ■ , A. We denote it as h\{%lj'\, \ip)\. Using 
Eq. (|T6l) . we may write the non- hermitian density as 

P0 + V5p{uj) = J2{\<P^){^^^\+V\M^)){^^\+V\<P^){y^iuJ)\} (27) 
i 

= Yl{\^^)+v\M^))}m\+v{Y^iu;)\}. (28) 

i 

In the last equation, we assume the linear approximation with respect to 1]. The fact that 



7 



Sh{uj) is proportional to Sp{uj) and Sh'^{uj) proportional to Sp^{uj) leads to 

ho + vShiu) = h[po + vSp{u;)] = h[{^\+r]{Y{u)\,\(P)+r]\X{u))], (29) 
ho + 7]6h\uj) = h[po + vSpH^)] = h[{<p\+v{X{uj)\,\<P)+v\Y{uj))], (30) 

up to the first order in a small parameter rj. In other words, the induced fields may be 
calculated using the finite difference with respect to r/; 

6hiuj) = ^{h[{^'l\^)]-hm,\<p)]), (31) 

where {ipi\ = + ri{Yi{u)\ and = + ri\Xi{u)). Its hermitian conjugate, 6h'^{uj), 
may be expressed as the same equation fl3T|) . but with {ip^\ = + ri{Xi{u)\ and = 
|0,) + r/|F,(^)). 

Using these numerical differentiation, the r.h.s. of the RPA equations, (12^ and (l25ll . 
can be easily calculated by action of the HF Hamiltonian, h[{tp'\, \-ip)~\ , on the single-particle 
orbitals, For the first sight, Eqs. and (l25l) do not look like linear equations. 

However, since 6h{uj) linearly depends on \Xi{uj)) and {Yi{uj)\, they are inhomogeneous linear 
equations with respect to \Xi{uj)) and {Yi{uj)\. It is obvious in a sense that they are equivalent 
to the matrix form of Eq. (flQl) . Therefore, we can employ a well-established iterative method 
for their solutions. If the linear equation is described by a hermitian matrix, the conjugate 
gradient method (CGM) is one of the most powerful method. However, in general, we may 
take the frequency uj complex, then, the RPA matrix becomes non-hermit ian. Then, we 
should use another kind of iterative solver, for instance, the bi-conjugate gradient method 
(Bi-CGM). A typical numerical procedure is as follows: (i) Fix the frequency uj that can 
be complex, and assume initial vectors {n = 0), \xl^\uj)) and {y}^\uj)\. (ii) Update the 
vectors, \xl"'^^\uj)) and {Y^"'^^\uj)\, using the algorithm of an iterative method, such as 
CGM and Bi-CGM. (iii) Calculate the residual of Eqs. ([21]) and ([25]). If its magnitude is 
smaller than a given accuracy, stop the iteration. Otherwise, go back to the step (ii) with 
n — > 77, + 1. 

The most advantageous feature of the present approach is that it only requires operations 
of the HF Hamiltonian, h[{ip'\, lip)]- These are usually included in computational programs 
of the static HF calculations. Only extra effort necessary is to estimate the HF Hamil- 
tonian with different bra and ket single-particle states, {ipH and Therefore, a minor 
modification of the static HF computer code will provide a numerical solution of the RPA 
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equations. Hereafter, we call this numerical approach "finite amplitude method". Appar- 
ently, the present method is also applicable to the RPA eigenvalue problems with a trivial 
modification. 



C. Transition strength in the linear response 

In this subsection, we present how to calculate transition strength using solutions of Eqs. 
(IMl) and (l25l) . Assuming that the system is at its ground state |$o) with energy Eq = 
at t = — oo, and that the external field V"ext(i^) is adiabatically switched on {u u ± ie in 
Eq. ([8])), the state at time t will be 

\^{t)) = |$o) -^5^e-^^"* r dtV^"*>„)(<l>„|Kxt(t')l*o). (32) 

in the first-order approximation with respect to Vcxt- Here, and En are the n-th excited 
state and its excitation energy, respectively. Especially, if the external field has a fixed 
frequency cu > 0, K;xt(^) = i]Fe~'^'^^ + r7*F^e*'^*, this is written as 

where F is an arbitrary one-body operator. Then, the expectation value of at time t is 

= (*o|f'l<l'o>+i)S(f;'^K"' + --- , (34) 

^ ' ^\iu-En + ie u + En-ie J ^ ^ 

Taking the limit of e ^ 0, we have the transition strength, 

^ V I {^n\F\%)mu - En) = --lmS{F; cu). (36) 

Ckj ^ — ^ TT 

n 

Comparing Eq. flMl) with the expectation value in the TDHF state, 

tr{FV(t)} = tr{FVo} + tr{F^5p(cu)}e"^^* + ■ • • . (37) 

S{F : a;) in the RPA is written as 

Sj^pa{F;uj) = tT{F^p{uj)} =tti{[Spy,Spiuj)]} (38) 
= J2 mFMM^)) + {Y^{uj)\F^\<l>^)) . (39) 

i 

Here, Spp is defined by Spp = i[F, po]. 



D. Separation of Nambu-Goldstone modes 

The RPA theory is known to have a property that the zero-energy modes are exactly 
decoupled from physical (intrinsic) modes of excitation. Since the zero modes are associated 
with the symmetry breaking in the HF ground state, it is also called "Nambu-Goldstone 
modes" (NG modes). When P is a hermitian symmetry operator of the Hamiltonian, 
[if, P] = 0, then, the transformed ground-state density, po = e*"^poe~*"^, also satisfies 
the HF equation (Q. Expanding the equation up to the first order in a, we have 

[ho, 5p] + [5h, Po] = 0, (40) 

where 

Sh 

Sp = Po- po = ia[P, Po], Sh = h[po] - ho = — ■ 5p. (41) 

dp 

This indicates that 6p is an RPA eigenmode corresponding to u = 0. 6p generated by the 
operator R conjugate to P {[R, P] = i) can be defined in a similar manner. For instance, the 
translational symmetry is expressed by the total momentum as P and the center-of-mass 
coordinate as R. We denote these transition densities associated with the NG mode as 

5pp ^ l[P,po] = -5p = Y, + \<i>^){P^\) , (42) 

i 

5pR = t[R,po] = J2 {\R^){<I>^\ + \<I^^){R^\) , (43) 
i 

where we have defined \Pi) = iP\(f)i) and \Ri) = iR\(l)i). Provided that P and R are 
hermitian, 5pp and 5pQ are also hermitian. Therefore, we cannot normalize them in terms 
of the normalization condition of Eq. ffTOj) . These modes are automatically orthogonal to 
other normal modes with 7^ 0. If we solve the RPA equation fully self-consistently, 
the NG modes should be clearly separated from other modes. However, in practice, we 
often encounter a mixture of spurious components in physical excitations. For instance, the 
coordinate-space is discretized in mesh to represent wave functions in Sec. Illlt which violates 
the exact translational and rotational symmetries. We also use a smoothing parameter F to 
make the frequency complex, then, low-lying excited states are embedded in a large tail of 
the NG-mode strength {5p{uj) 00 for ^ ). Here, we present a prescription to remove 
the strength associated with the NG mode. 

Let us assume that there is a mixture of NG modes in a calculated transition density. 
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5pcai(t^) = 5pphy(t^) + XpSpp + Xr6pr, (44) 

where "physical" transition density, 5pphy(ci;), is free from the NG modes. Here, we assume 
there is a single NG mode, for simplicity. It is straightforward to extend the present pre- 
scription to the one for more than one NG modes. Since (5pphy should be orthogonal to the 
NG modes, (5pphy should satisfy 

{%\[6pp,6pp^y{uj)]\%) = {%\[6pR,6ppi,y{uj)]\<l>o) = 0. (45) 

Utilizing the canonicity condition, [R,P] = i, the orthogonality condition, Eq. fH51) . deter- 
mines the coefficients, Xp{R), as 

Ap = - {Y^i^)\R^)) , (46) 

i 

i 

Substituting these into Eq. (jBl), we may extract SppYiy{uj) from the "contaminated" transi- 
tion density Spcaii^^)- 

III. NUMERICAL APPLICATIONS 
A. Coordinate-space representation 

In case of zero-range effective interactions, such as Skyrme interactions, the HP Hamil- 
tonian, h{r) = h[p{r)], is a functional of local one-body densities. Then, it is convenient 
to adopt the coordinate-space representation. In the foUowings, we assume r involves the 
spin and isospin indexes, if necessary. The RPA equations, and (|25|) . for a complex 
frequency u can be written in the r-representation as 

{ho{r) - ei - uj)X,{r,uj) + Sh{r,uj)(j)i{r) = -1/ext(r, cj)0i(r), (48) 
{{ho{r)-e, + u;*)Y,{r,u;) + 6h\r,u;)<Pi{r)y = - { KL(r, u;)0i(r)}* . (49) 

Here, for simplicity, we omit the projection operator, Q, on both sides of these equations. 
In the finite amplitude method, the operation of 6h{r, u) is calculated by 

5h{v,cu)^,{r) = ^ {h[^'*,^]iv)<P,iv) - e,0,(r)) , (50) 
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with 4''*{r) = 0*(r) + riY*{uj, r) and V'j(r) = + '7-^i(i"7 Exchanging the forward and 
backward amphtudes in ipi{r) and ipi{r), we may calculate Sh'^ {r,uj)(f)i{r) in the same way. 
Adopting the fixed-a; local external field 

V,,,iT,u')=S^^,FiT), (51) 

the transition strength can be obtained from the calculated forward and backward ampli- 
tudes, 

^ J2\{n\F\0)mu;-E^), (52) 

n 

= -l^'^T. j dv{<P:{v)F\v)X,{v,uj) + Y:{v,u)F\v)<P,{v)]. (53) 

i 

We apply the present method to the BKN interaction which contains two-body (zero- and 
finite-range) and three-body interactions. For this schematic interaction, the spin-isospin 
degeneracy is assumed all the time and the Coulomb potential acts on all orbitals with a 



charge e/2 13|. The HF one-body Hamiltonian in the coordinate-space representation is 
given by 

h[p] = + ^toP(r) + ^p^(r) + Wy[p]{v) + l^c[p](r), (54) 

where the Yukawa potential, Wy, and Coulomb potential, Wc, consist of their direct terms 
only. For the finite amplitude approach, it is convenient to rewrite Eq. as 

. . A/4 r A/4 

i=l I 1=1 

I- A/A 

+ / rfr't;(r-r')^4^i(r')C(r'), 

i=i 



(55) 



where f (r) is a sum of the Yukawa and the Coulomb potential. 



^(r) ^ V,a- + (56) 

r r 



We adopt the parameter values from Ref. 



B. Numerical details 



We use the three-dimensional (3D) coordinate-space representation for solving the RPA 
equations. The model space is a sphere of radius of 10 fm, discretized in square mesh 
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of Ax = Ay = Az = 0.8 fm. The number of grid points in the sphere is 8217. The 
differentiation is approximated by the nine-point formula. The frequency u is varied from 
zero to 40 MeV with a spacing of Au = 200 keV (201 points). A small imaginary part is 
added to u: u uj + ir/2 with F = 500 keV. In numerical calculations, we use real variables 
with double precision (8 bytes) and complex variables of 8 x 2 bytes. In Eq. flSOp . we choose 
the parameter rj in '?/'i(r) = 0i(r) + riXi{r) and '?/'j'*(r) = 0*(r) + riY*{r), as follows: 

1Q-5 1 I 

"- ma.iiVm.TVy)} - ^ (67) 

In order to obtain the forward and backward amplitudes at a frequency uj, we adopt the 
Bi-CGM as an iterative solver for Eqs. fj48l) and fH9|) . starting from the initial values of 
Xj(r) = Y*{r) = 0. We set the convergence condition that the ratio of the remaining 
difference to the r.h.s. of Eqs. ( HHl) and ( l49l) is less than 10^^. The number of iteration 
necessary to reach the convergence depends on the choice of the external field Vext{^), 
the frequency u, the smoothing parameter F, and residual interactions included in the 
calculation. The convergence is relatively slow for an external field coupled to the NG modes. 
A larger number of iteration is required for a larger u value. Typically, the calculation 
reaches the convergence in 10 to 100 iterations for uj < 10 MeV, but it requires more than 
500 iterations for uj > 30 MeV. The number also depends on the smoothing parameter 
F. Roughly speaking, larger number of iteration seems to be required for smaller F. If 
we neglect the residual Coulomb and Yukawa interactions of finite range, the convergence 
becomes much faster. We solve the differential equations to obtain the Coulomb and Yukawa 



potentials using the CGM jl5 |. 

VVc = -27reV(r), (^\/^ - 1"^ Vy = -i7iVoap{r) . (58) 

It turns out to be important to solve these equations with high accuracy. We set the 
convergence condition that the ratio of the remaining difference to the r.h.s. of Eq. (1581) is 
less than 10"^'^. Since the convergence of the CGM is very fast, this is not a problem. 

C. Results 

In this section, we show calculated response for isoscalar (IS) modes of compressional 
dipole, quadrupole and octupole for ^°Ne. The main purpose of the calculation is to test 
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capability of the present numerical approach, the finite amplitude method. The ^"Ne nu- 
cleus has a prolate shape with a quadrupole deformation (3 ~ 0.4 in the HF ground state. 
Identifying the symmetry axis with z-axis, we use external fields with a fixed frequency, 
Vext(r,cj) = (5Ax(r), 



r2r2o(r), r^Y^iii), r^Y22{r), 
^r=^r3o(r), r3r3i(r), r^Y:i2{i), 
Then, the strength distribution, 

dB{uj\QxK) 



for IS dipole, A = 1 

for IS quadrupole, A = 2 (59) 



r^y33(f), for IS octupole, A = 3. 



(60) 



will be calculated according to Eq. fl53l) . 



1. Isoscalar quadrupole response: Accuracy of the finite amplitude method 

In Fig. [T](a), we show results for the IS quadrupole strength distribution. There is a 
NG mode in the K = 1 sector, corresponding to the nuclear rotation. This is clearly 
seen in the response of the K = 1 mode, having a large peak near uj = 0. The RPA 
correlation brings the lowest one-particle-one-hole (Iplh) excitation at = 4.5 MeV down 
to zero. The response function for the K = 1 mode was not obtained by the small-amplitude 
TDHF method in Ref. because the nucleus actually rotates in real time, that violates 
the small-amplitude approximation. This is an advantage of the present method over the 
time-dependent approach. The lowest intrinsic (physical) excitation corresponds to the 
K = 2 mode at u = 8 MeV which is close to energy of the Iplh excitation. This suggests 
that the correlation effect is weak for this mode, supported by a small K = 2 quadrupole 
strength at = 8 MeV. In contrast, the next lowest mode at Ex = 9.6 MeV with = is 
somewhat lowered by the correlation and exhibit a larger strength. Ref. |l^ shows results 
of configuration mixing calculation with the BKN interaction, indicating = 0+ around 
Ex = 7 MeV and J"" = 2+ state near 8 MeV. Large peaks at = 15 ~ 22 MeV should 
correspond to the IS giant quadrupole resonance. It clearly shows deformation splitting; the 
K = peak at lowest, the K = 1 in the middle, and the K = 2 at the highest energy. 
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FIG. 1: IS quadrupole strength distribution for ^'^Ne. The sohd, dotted, dashed hnes indicate those 
with K = 0, 1, and 2, respectively. Results of two kinds of calculations, (a) the finite amplitude 
method (FAM) and (b) the conventional RPA, are compared. 

Now, let us demonstrate accuracy of the present finite ampUtude method. In Fig. [H 
results of the conventional RPA, which explicitly estimates the residual interactions 6h/Sp, 
are also presented in the panel (b). These two kinds of calculations, (a) and (b), provide 
identical results in the accuracy of three to four digits. 



2. Isoscalar dipole and octupole response: Removal of NG modes 



Next, we show the strength distribution for the isoscalar compressional dipole mode. This 
mode has been of significant interest because its energy is related to the compressibility 
of nuclear matter, providing information independent from the monopole resonance. The 
compressional modes in spherical nuclei have been extensively studied with the continuum 



RPA calculations 
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19|. However, these calculations are not fully self-consistent, thus, 



need to remove mixture of the NG (translational) components by modifying the dipole 
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FIG. 2: IS compressional dipole strength distribution for ^'^Ne. Strengths associated with the K = 
modes are shown in the left (a) and those with K = 1 in the right (b). The upper panels show 
strengths calculated with 6pcaii^), while the lower panels, (a-3) and (b-3), show those calculated 
with 5pp\iy{uj) in Eq. (jH]). See text for details. 



operator. This produces some ambiguity in their results. In fact, the importance of the 



2l|. So far, our 



full self-consistency has been stressed for the compressional modes [20|, 
understanding of the compressional dipole mode is still obscure and further studies are 
needed. In this section, we show a fully self-consistent calculations for deformed nuclei. 
In Fig. [21 the compressional dipole strength is shown, K = mode at the left (a) and 
= 1 at the right (b). There are the NG modes associated with the translational symmetry 
breaking near uj = 0, seen in Figs. [2] (a-1) and (b-1). These NG peaks are so huge that other 
peaks are invisible in these insets. The vertical axis is magnified in the panels (a-2) and 
(b-2). The giant resonance peaks are spread over = 16 ~ 30 MeV for K = and 20 — 40 
MeV for K = 1. There is a sharp peak at uj = 4.5 MeV, which is embedded in the tail of 
the NG mode. In order to estimate the strength carried by this state, we need to separate 
out the contribution from the NG mode. This is done by using the prescription described 
in Sec. Ill D\ adopting the center-of-mass coordinates and the total linear momenta (3 NG 
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FIG. 3: IS octupole strength distribution for ^''Ne. Strengths associated with the K = 0, 1, 2, and 
3 octupole modes are shown in the panels (a), (b), (c), and (d), respectively. In the panels (a) and 
(b), the strengths calculated with 5pca.i{uj) are presented by dotted lines, while those with 6pphy{^) 
are by solid lines. 

modes). Strength associated with the "physical" transition density (5pphy(to') is shown in 
Fig. [2] (a-3) and (b-3). The large strength of translational modes near = is properly 
removed. The other physical peaks with finite u are unchanged, which indicates that there 
is very little mixture of the NG modes because our calculation is fully self-consistent. Now, 
we may identify the K = peak at = 4.5 MeV as an isolated peak. 

Finally, we show IS octupole strength distribution with K = 0, 1, 2, and 3 in Fig. [31 The 
lowest octupole state is at = 4.5 MeV with K = 0, and the second lowest is aX u = 8.1 
MeV with K = 2. These results are similar to that of the variation-after-parity-projection 
calculation 22 1 and that of the configuration- mixing calculation [l^. Experimentally, the 
band head of the K = 2 band {J'^ = 2^) is observed at = 5.0 MeV and that of i^' = 
{J^ = 1~) is at = 5.8 MeV. The BKN interaction, that does not contain the spin-orbit 
force, is able to reproduce the K = state in a reasonable accuracy, however, fails to provide 
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a quantitative description for the K = 2 state. This suggests that the spin-orbit force does 
not play an important role for the K = state. In fact, the parity-projected HF calculation 
with the Skyrme interaction has confirmed very small contribution of the spin-orbit force in 
this K"" = 0~ state [23]. 

Since the nucleus is deformed, the dipole modes are coupled to the octupole modes. We 
may identify peaks at the same positions in Figs. [2] and [3] for the K = and K = 1. We see 
a small spurious K = 1 peak near uj = 0, even after removing the NG components (solid line 
in Fig. [3](b)). However, the peak height of the NG mode is about 5,000 fm^/MeV. Thus, 
more than 98 % of the NG strength is actually removed. We can say that the method in 
Sec. IIIDI also works for octupole modes. 

IV. CONCLUSION AND DISCUSSION 

We have presented a new numerical approach to the RPA calculation, "finite ampli- 
tude method". The finite amplitude method does not require complicated programming for 
complex residual interactions. Instead, it resorts to the numerical derivation of the residual 
interaction (induced field), 6h/6p-6p. The most advantageous feature of the present method 
is its feasibility of programming a computer code. The RPA calculation can be accomplished 
with a minor extension of the static HF computer code, to construct the HF Hamiltonian 
with independent bra and ket single-particle states. 

Here, we would like to make a remark on the meaning of different bra and ket states. 
This does not mean matrix elements between different Slater determinants which are rather 
complicated. These "off-diagonal" elements are necessary for configuration-mixing calcu- 
lations, such as the generator-coordinate method. The finite amplitude method does not 
require these. All we need is "diagonal" matrix elements of a certain one-body operator, 6, 
in the linear order with respect to variation of the single-particle states, 

AAA A 

J2{i^^\0\A) ~ Y.^<P^\^\<P,) + J2{^^^\0\<P^) + (61) 
i=l 1=1 1=1 1=1 

In order to calculate the second and the third terms in the r.h.s. separately, we need to input 
independent bra and ket single-particle states. This can be achieved by a minor extension 
of the static HF code. 

The method has been applied to calculations of the isoscalar dipole, quadrupole, and 
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octupole response functions. Since the adopted interaction is rather schematic, we do not 
discuss here calculated properties of these modes. Instead, we would like to emphasize 
characteristic features of the finite amplitude method. First of all, the transition density 
coupled to the NG modes can be calculated without any special treatment. As is seen in 
the quadrupole {K = 1) and octupole modes {K = and 1) in Figs. [T] and [31 the NG modes 
appear in responses to a variety of external fields, especially for deformed nuclei. This causes 
a serious problem in the time-dependent calculation of the small-amplitude TDHF [ll^. If 
we want to remove the strength associated with the NG modes, we can use the prescription 
presented in Sec. IIIDI This has been demonstrated in the dipole and octupole strength 
distributions. Second, since we do not calculate the residual interaction explicitly, it is easy 
to carry out the fully self-consistent RPA calculation for realistic interactions including spin- 
orbit, derivative, and Coulomb terms. The implementation of the present method does not 
depend on the complexity of the interactions. For instance, the compressional dipole mode 
has been a long-standing problem in microscopic calculations 2J]. The problem is related 
to difficulties in the fully self-consistent treatment and the coupling to the translational 
modes. Our new approach may provide a tool to clarify this point. Last but not least, the 
finite amplitude method is an efficient method to calculate the strength distribution. We 
may control the necessary energy resolution by the smoothing parameter F. The numerical 
application to the BKN functional shows that its efficiency is next to the time- dependent 
method, better than the other methods including the Green's function method HI] and the 
diagonalization method jl4|. The diagonalization of the RPA matrix is very efficient if we 
are interested in only a few lowest states, however, it becomes more and more difficult for 
higher excitation energies. 

For future developments, it is interesting to combine the present method with the 
absorbing-boundary-condition approach in Ref. ll|]. This enables us to calculate response 
in the continuum, overcoming difficulties in the time-dependent method. It is also very in- 
teresting to extend the method in the HF scheme to the one in the Hartree-Fock-Bogoliubov 
framework. In this paper, we adopt a simple interaction to check the method, but the finite 
amplitude method shows its real power for a complex density functional. Applications of 
the method to the realistic Skyrme functionals are under progress at present and will be 
published in near future. 
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